De l’assurance qualité à l’assignation taxonomique
2025-10-22
Note
La présentation est disponible à l’adresse: https://taq.info/trainingMetabarcoding/
Note
Le signal renvoyé par le laser est traduit en chromatographe pour chaque amplicon séquencé.
Source: LabXchange.org
Important
Le score est basé sur:
Les sortie du séquenceur Illumina
Exemple de nomenclature de fichier de sortie du séquenceur Illumina: sample1_1_L001_R1_001.fastq.gz
| Signification | Exemple | |
|---|---|---|
| sample1 | Identifiant de l’échantillon | sample1, sample2 |
| _1 | Numéro de réplicat | 1, 2, 3… |
| L001 | Numéro de ligne / lane: | L001-L008 |
| R1/R2 | Direction de lecture | R1=avant, R2=arrière |
| 001 | Segment de fichier | 001, 002… |
Note
Contient plusieurs lignes dont un ensemble de 4 lignes correspond à un amplicon.
system.file(package="ShortRead", "extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz") |>
gzfile() |>
readLines(n = 12) [1] "@ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1"
[2] "GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGACATGAAGTCAATGAAGGCCTGGAATGTCACTACCCCCAG"
[3] "+"
[4] "HHHHHHHHHHHHHHHHHHHHEBDBB?B:BBGG<DDAA?AABFEFBDBD@DDECEE3>:?;@@@>?=BAB?##"
[5] "@ERR127302.21406531 HWI-EAS350_0441:1:88:9330:2587#0/1"
[6] "CTAGGGCAATCTTTGCAGCAATGAATGCCAATGGGTAGCCAGTGGCTTTTGAGGCCAGAGCAGACCTTCGGG"
[7] "+"
[8] "IIIIHIIIGIIIIIIIHIIIIEGBGHIIIIHGIIHIIIIIIIHIIIHIIIIIGIIIEGIIGBGE@DDGGGIG"
[9] "@ERR127302.22173106 HWI-EAS350_0441:1:91:10434:14757#0/1"
[10] "TGGGCTGTTCCTTCTCACTGTGGCCTGACTAAAACCCAGTGGCATTAAGAAAGAGTCACGTTTCCCAAGTCT"
[11] "+"
[12] "GGHBHGBGGGHHHHDHHHHHHHHHFGHHHHHHHHHHHHHHHHHHHHHGHFHHHHHHHHHHHHHH8AGDGGG>"
system.file(package="ShortRead", "extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz") |>
gzfile() |>
readLines(n = 4)[1] "@ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1"
[2] "GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGACATGAAGTCAATGAAGGCCTGGAATGTCACTACCCCCAG"
[3] "+"
[4] "HHHHHHHHHHHHHHHHHHHHEBDBB?B:BBGG<DDAA?AABFEFBDBD@DDECEE3>:?;@@@>?=BAB?##"
@)+)Décomposition de l’identifiant
@SH00321:6:BWR98207-2813:1:1101:1065:1015 1:N:0:AACCATAGAA+GGCGAGATGG
Exemple de chaine de caractères donnant le score de qualité (1 caractère par nucléotide): GGGGGGGGGGGGGGGGGGGGGGGG9GG-G9G9G9GGGG
| Caractère ASCII | Score Phred | Taux d’erreur | Précision |
|---|---|---|---|
G |
38 | 0,016% | 99,984% |
9 |
24 | 0,40% | 99,60% |
- |
12 | 6,31% | 93,69% |
Tip
Règle générale : Q30 ou plus = haute qualité (99,9% de précision)
Les scores de qualité utilisent l’encodage ASCII.
Quality encoding: !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHI
| | | | |
Quality score: 0........10........20........30........40
Avec R, comment convertir l’encodage ASCII en score de Phred et en probabilité d’obtenir une erreur (précision)?
# Convertir la chaîne de qualité en scores Phred
quality_string <- "GGGGGGG99GG-GGG"
quality_scores <- as.integer(charToRaw(quality_string)) - 33
# Afficher la correspondance
data.frame(
Caractere = strsplit(quality_string, "")[[1]],
Score_Phred = quality_scores,
Precision = paste0(round((10^(-quality_scores/10)) * 100, 2), "%")
) Caractere Score_Phred Precision
1 G 38 0.02%
2 G 38 0.02%
3 G 38 0.02%
4 G 38 0.02%
5 G 38 0.02%
6 G 38 0.02%
7 G 38 0.02%
8 9 24 0.4%
9 9 24 0.4%
10 G 38 0.02%
11 G 38 0.02%
12 - 12 6.31%
13 G 38 0.02%
14 G 38 0.02%
15 G 38 0.02%
[1] 20000
# Extraire les amplicons
sequences <- ShortRead::sread(fq)
# Extraire les identifiants
ids <- ShortRead::id(fq)
# Extraire les scores de qualité
qualities <- Biostrings::quality(fq)
# Premier identifiant et séquence
ids[1]BStringSet object of length 1:
width seq
[1] 53 ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1
DNAStringSet object of length 1:
width seq
[1] 72 GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGAC...GTCAATGAAGGCCTGGAATGTCACTACCCCCAG
Min. 1st Qu. Median Mean 3rd Qu. Max.
72 72 72 72 72 72
# Distribution des scores de qualité
qual_matrix <- as(qualities, "matrix")
mean_quality <- rowMeans(qual_matrix)
summary(mean_quality) Min. 1st Qu. Median Mean 3rd Qu. Max.
2.00 33.53 37.35 34.84 38.85 39.99
A C G T
[1,] 0.1944444 0.2638889 0.2916667 0.2500000
[2,] 0.2361111 0.2222222 0.3194444 0.2222222
[3,] 0.2361111 0.2638889 0.2222222 0.2777778
[4,] 0.2500000 0.2638889 0.1666667 0.3194444
RqcNote
Exercice simple: Utilisez le code R que nous venons de voir pour analyser vos propres fichiers FASTQ
data/RqcTip
Si vous n’avez pas de données, vous pouvez télécharger l’exemple ici. Ce sont 10 échantillons de métabarcoding mitofish-12S
10:00
Hakimzadeh et al, 2024
cut, sort, uniq etc. qui manipulent des chaînes de caractères ou des nombresVue d’ensemble du pipeline de métabarcoding
SampleID_*_R1_001.fastq.gzPour chaque lecture fusionnée:
Télécharger BARQUE en cliquant sur le lien
Ou en rendant à l’adresse suivante: https://github.com/enormandeau/barque
fastq.gz avec la bonne nomenclature dans le dossier 04_data02_info/primers.csv02_info/barque_config.shTip
Si vous voulez réutiliser la configuration du pipeline pour d’autres fois, c’est une bonne pratique de dupliquer le fichier et de le sauvegarder sous un autre nom.
04_dataOn peut faire tout de suite l’étape 1 avec vos données.
Tip
Si vous avez pas de données, vous pouvez télécharger celles fournies par Étienne Normandeau, en cliquant ici. Ce sont 10 échantillons de métabarcoding mitofish-12S, chacun.
02_info/primers.csv# pour la ligne correspondante à votre amorce# sur les autres lignes# Modify the following parameter values according to your experiment
# Do not modify the parameter names or remove parameters
# Do not add spaces around the equal (=) sign
# It is a good idea to try to run Barque with different parameters
# Global parameters
NCPUS=20 # Number of CPUs to use. A lot of the steps are parallelized (int, 1+)
PRIMER_FILE="02_info/primers.csv" # File with PCR primers information
# Skip data preparation and rerun only from vsearchp
SKIP_DATA_PREP=0 # 1 to skip data preparation steps, 0 to run full pipeline (recommended)
# Filtering with Trimmomatic
CROP_LENGTH=200 # Cut reads to this length after filtering. Just under amplicon length
# Merging reads with flash
MIN_OVERLAP=30 # Minimum number of overlapping nucleotides to merge reads (int, 1+)
MAX_OVERLAP=280 # Maximum number of overlapping nucleotides to merge reads (int, 1+)
# Extracting barcodes
MAX_PRIMER_DIFF=8 # Maximum number of differences allowed between primer and sequence (int, 0+)
# Running or skipping chimera detection
SKIP_CHIMERA_DETECTION=0 # 0 to search for chimeras (RECOMMENDED), 1 to skip chimera detection
# or use already created chimera cleaned files
# vsearch
MAX_ACCEPTS=20 # Accept at most this number of sequences before stoping search (int, 1+)
MAX_REJECTS=20 # Reject at most this number of sequences before stoping search (int, 1+)
QUERY_COV=0.6 # At least that proportion of the sequence must match the database (float, 0-1)
MIN_HIT_LENGTH=100 # Minimum vsearch hit length to keep in results (int, 1+)
# Filters
MIN_HITS_SAMPLE=10 # Minimum number of hits in at least one sample to keep in results (int, 1+)
MIN_HITS_EXPERIMENT=20 # Minimum number of hits in whole experiment to keep in results (int, 1+)
# Non-annotated reads
NUM_NON_ANNOTATED_SEQ=200 # Number of unique most-frequent non-annotated reads to keep (int, 1+)
# Multiple hits
MIN_DEPTH_MULTI=10 # Min depth to report unique reads per sample in multiple hit reports
# OTUs
SKIP_OTUS=1 # 1 to skip OTU creation, 0 to use it
MIN_SIZE_FOR_OTU=20 # Only unique reads with at least this coverage will be used for OTUsOn ajuste les paramètres du pipeline dans le fichier 02_info/barque_config.sh
Tip
Pour déterminer le nombre de coeur sur votre ordinateur, utiliser la commande R suivante parallel::detectCores() - 2
BARQUE utilise différentes bases de données selon le type d’amorces utilisé pour l’amplification:
Bases de données intégrées:
Formats requis:
03_databases/Les bases de données BOLD et SILVA pré-formatées pour BARQUE sont disponibles ici:
https://www.ibis.ulaval.ca/services/bioinformatique/barque_databases/
Important
Après téléchargement, renommer le fichier en silva.fasta.gz ou ****.fasta.gz** et le placer dans 03_databases/
VSEARCH utilise une approche de similarité de séquences basée sur l’alignement global:
Dans le fichier barque_config.sh, ces paramètres contrôlent l’assignement:
MAX_ACCEPTS=20 # Nombre maximal de correspondances acceptées avant d'arrêter
MAX_REJECTS=20 # Nombre maximal de rejets avant d'arrêter la recherche
QUERY_COV=0.6 # Couverture minimale de la requête (60% de la séquence)
MIN_HIT_LENGTH=100 # Longueur minimale d'alignement pour conserver le résultatTip
QUERY_COV=0.6 signifie qu’au moins 60% de votre séquence doit correspondre à une séquence de la base de données pour être considérée comme un match valide
Séquence inconnue:
ATCGATCGATCGATCG...
(150 nucléotides)
Meilleure correspondance dans la base:
Critères de validation:
Résultat: L’amplicon est assigné à Salmo salar avec un niveau de confiance élevé
Que se passe-t-il si plusieurs espèces ont le même score?
Warning
Les correspondances multiples peuvent indiquer:
BARQUE requière plusieurs programmes qui ne sont pas tous compatibles avec l’environnement Windows
Hakimzadeh et al, 2024
Docker crée un environnement isolé (conteneur) qui contient:
Le conteneur peut accéder aux fichiers de votre ordinateur via des volumes montés.
Tip
Docker permet d’utiliser BARQUE sur Windows, Mac ou Linux sans installer manuellement tous les programmes requis
@Soroosh Nazem
# Base image with Ubuntu
FROM ubuntu:24.04
# Update and install essential packages
RUN apt-get update && apt-get install -y \
wget \
curl \
build-essential \
software-properties-common \
default-jre \
parallel \
r-base \
bc \
python3-setuptools \
ca-certificates \
&& rm -rf /var/lib/apt/lists/*
# Install FLASH (v1.2.11+)
RUN wget https://github.com/dstreett/FLASH2/archive/refs/tags/2.2.00.tar.gz && \
tar -xzf 2.2.00.tar.gz && \
cd FLASH2-2.2.00 && \
make && \
mv flash2 flash && \
cp flash /usr/local/bin && \
cd .. && rm -rf FLASH2-2.2.00 2.2.00.tar.gz
# Install VSEARCH v2.14.2+
RUN wget https://github.com/torognes/vsearch/releases/download/v2.30.0/vsearch-2.30.0-linux-x86_64.tar.gz && \
tar -xzf vsearch-2.30.0-linux-x86_64.tar.gz && \
cp vsearch-2.30.0-linux-x86_64/bin/vsearch /usr/local/bin && \
rm -rf vsearch-2.30.0-linux-x86_64*
# Default shell
CMD ["/bin/bash"]Ouvrir un terminal (MacOS) ou PowerShell (Windows) exécuter la commande suivante pour télécharger l’image BARQUE:
Cette commande télécharge l’image Docker contenant BARQUE et tous ses programmes requis (FLASH, VSEARCH, Trimmomatic, etc.)
Note
Le téléchargement peut prendre quelques minutes selon votre connexion Internet
Pour lancer BARQUE dans un conteneur Docker avec accès à vos fichiers:
# Pour macOS
docker run -it -v /chemin/vers/barque:/barque steveviss/barque:latest
# Pour Windows
docker run -it -v C:\chemin\vers\barque:/barque steveviss/barque:latestImportant
Remplacer /chemin/vers/barque par le chemin absolu vers votre dossier BARQUE local
Explication des options:
-it : Mode interactif avec terminal-v /chemin/vers/barque:/barque : Monte votre dossier local dans le conteneursteveviss/barque:latest : L’image Docker à utiliserUne fois dans le conteneur Docker, vous êtes dans un environnement Linux avec BARQUE prêt à l’emploi:
Application disponible à l’adresse suivante: https://cloud.taq.info, contacter moi pour un accès: steve.vissault@inrs.ca
Application disponible à l’adresse suivante: https://cloud.taq.info, contacter moi pour un accès: steve.vissault@inrs.ca
Vous pouvez installer la version de développement de barqueReport depuis GitHub avec :
La fonction principale génère un rapport HTML interactif à partir des résultats du pipeline BARQUE:
library(barqueReport)
# Sélection interactive du dossier BARQUE et du fichier de qualité Illumina
generate_metabarcoding_report(
barque_output_folder = file.choose(), # Sélectionner le dossier de sortie barque
samples_ids = c("ST1", "ST2", "ST3", "ST4"),
title = "Étude du Lac Supérieur",
subtitle = "",
author = "Steve Vissault, Julie Couillard & Tuan Ahn To",
illumina_quality_file = file.choose(), # Sélectionner le fichier Quality_Metrics.csv
blank_lab = "LAB_BLANK",
blank_field = "FIELD_BLANK"
)
Comment le séquencage est effectué?